In the oil and gas industry, seismic prospecting techniques are commonly used to aid in the search for and evaluation of subterranean hydrocarbon deposits. In seismic prospecting, a seismic source is used to generate a seismic signal that propagates into the earth and is at least partially reflected by subsurface seismic reflectors (i.e., interfaces between underground formations having different elastic properties). The reflected signals (known as "seismic reflections") are detected and recorded by seismic receivers located at or near the surface of the earth, in an overlying body of water, or at known depths in boreholes, and the resulting seismic data may be processed to yield information relating to the subsurface formations.
Seismic prospecting consists of three separate stages: data acquisition, data processing, and data interpretation. The success of a seismic prospecting operation depends on satisfactory completion of all three stages.
The seismic energy recorded by each seismic receiver during the data acquisition stage is known as a "seismic data trace". Seismic data traces typically contain both the desired seismic reflections and one or more unwanted noise components which can obscure or overwhelm the seismic reflections. One of the primary objectives of the data processing stage is to remove or at least attenuate these unwanted noise components so that the desired seismic reflections can be clearly identified and interpreted.
One method for attenuating unwanted noise components in seismic data traces is through the common-midpoint (CMP) stacking process. As will be well known to persons skilled in the art, the "midpoint" for a seismic data trace is the point midway between the source location and the receiver location for that trace. According to the CMP method, the recorded seismic data traces are sorted into common-midpoint gathers each of which contains a number of different seismic data traces having the same midpoint but different source-to-receiver offset distances. The seismic data traces within each CMP gather are corrected for statics (i.e., the effects of variations in elevation, weathered layer thickness and/or velocity, or reference datum) and normal moveout (i.e., the variation of traveltime with respect to source-to-receiver offset) and are then summed or "stacked" to yield a stacked data trace which is a composite of the individual seismic data traces in the CMP gather. Typically, the stacked data trace has a significantly improved signal-to-noise ratio compared to that of the unstacked seismic data traces.
As is well known in the art, for a horizontally layered earth having a single horizontal reflector, seismic signal traveltime bears a hyperbolic relation to source-to-receiver offset distance, as follows: ##EQU1## where "t" is the seismic signal traveltime for source-to-receiver offset "x"; "t.sub.0 " is the seismic signal traveltime for zero offset; and "V.sub.NMO " is the velocity of the seismic signal (commonly referred to as the "normal-moveout velocity" or "NMO velocity").
In the CMP stacking process, seismic data traces from different source-to-receiver offsets within a CMP gather are stacked by moving them along hyperbolas defined by equation (1) to the zero-offset position and then adding them together. As is well known in the art, the NMO velocity for a CMP gather is not a constant. Typically, the NMO velocity increases sporadically as two-way, zero-offset traveltime increases. For this reason, proper stacking of the seismic data traces within a CMP gather requires knowledge of the NMO velocity as a function of two-way, zero-offset traveltime.
Typically, determination of the NMO velocity function is done manually by expert seismic analysts. One conventional method for determining the NMO velocity function is through the use of velocity spectra. According to this method, the individual seismic data traces in a CMP gather are repeatedly NMO-corrected and stacked using a range of trial velocity values. The resulting stacked traces are then displayed side-by-side on a plane of velocity versus two-way, zero-offset traveltime, known as the "velocity spectrum". The velocity which results in the highest stacked amplitude for a given reflection is selected or "picked" as the NMO velocity for that reflection. The NMO velocity function may then be expressed as a set of velocity-traveltime pairs.
Many other conventional methods for determining or "picking" the NMO velocity function are known. Typically, these methods are based on determining the velocity that corresponds to the best coherency of the signal along a hyperbolic trajectory over the entire CMP gather. Coherency may be measured in a variety of ways, such as an unnormalized crosscorrelation, a normalized crosscorrelation, and energy-normalized crosscorrelation, or semblance (i.e., stack energy normalized by the mean energy of the individual traces going into the stack).
Further information on the use of velocity spectra and other known methods of velocity analysis may be found in Yilmaz, O., Seismic Data Processing, 1987, Society of Exploration Geophysists, Tulsa, Okla., pp. 166-183.
The accuracy of the NMO velocity functions resulting from these prior art methods is limited by the discrete step size used in the set of trial velocity values. As would be well known to persons skilled in the art, manual velocity picking by even the most experienced analyst can only be accurate to within a few percent, and even a small error in the NMO velocity can result in a significant loss in signal-to-noise ratio and bandwidth in the stacked data. Furthermore, the phenomenon of "stretching" (i.e., the change in wavelet shape produced by applying a normal-moveout correction) further degrades the measurement of coherency for the coherency-based methods. Muting can be used to minimize the effects of stretching, but this reduces the fold of the stacking process (i.e., the number of individual traces that contribute to the final stacked trace for the CMP location).
In conventional seismic data processing, the NMO velocity function is determined by a seismic analyst at a plurality of laterally spaced-apart locations along the seismic line. This is a time consuming process, and, for practical reasons, the number of locations where the NMO velocity function determination is made must be limited. Stacking programs then interpolate these analyst-determined NMO velocity functions to determine the NMO velocity function for each intermediate CMP location.
Even if the analyst-determined NMO velocity functions are correct, it is not certain that the process of interpolation gives the correct NMO velocity function at intermediate CMP locations. Typically, linear interpolation is used. Lateral variations in seismic velocity need not follow a linear variation. Even small errors (e.g., one to two percent) in the NMO velocity function used for stacking can degrade the quality of the stack. This problem would be familiar to persons skilled in the art.
Obviously, it is desirable to have an accurate NMO velocity function for every CMP location along a seismic line. Moreover, given the vast amount of data resulting from a typical three-dimensional (3-D) seismic survey, a method for automatically determining accurate NMO velocity functions is highly desirable.
Prior art attempts to develop automatic methods for determining the NMO velocity function have met with only limited success. In a typical CMP gather, there may be multiple reflections and refractions in addition to the primary reflections of interest. These multiple reflections and refractions may produce spurious peaks in the coherence spectra making it difficult to determine the peak that corresponds to the correct NMO velocity. An experienced seismic analyst often can disregard these spurious peaks and correctly pick the NMO velocity function. However, a completely automatic velocity picking procedure sometimes produces an unsatisfactory solution because the multiple reflections or refractions may have greater energy than the primary reflections of interest.
Some prior art automatic velocity picking procedures attempt to emulate the main criteria followed by seismic analysts by introducing constraints on the allowable interval or stacking velocities. Starting with an initial estimate of the velocity function, a conjugate gradient scheme is used to give an improved estimate of the velocity spectra. As a starting point, this procedure requires computation of the coherency over a range of zero-offset times and velocities for a number of different CMP locations. The choice of the initial estimate can, under certain conditions, lead to different results, and the effects of stretching are present. In any case, the resolution of these automated procedures is no finer than the step size used in the coherency computation.
The prior art also teaches the addition of a penalty function to the conventional coherency measurement; this function penalizes velocity models that differ in smoothness from an initial interval-velocity model. This requires computation of the coherency over a range of zero-offset times and velocities for a number of different CMP locations. The resolution of this method is limited to the interval size used in defining the trial velocities and times, and stretching remains a problem.
Besides the limit in resolution due to the size of the velocity or time step used in the computation, the coherency-based methods are also subject to an additional loss of resolution because the coherency between the traces in a CMP gather is computed over a time window which introduces a smearing of the velocity function. This smearing can be particularly problematic when the window encompasses a strong reflection on the seismic trace. Shifting the center of the window will often produce no discernible change in the coherency value. Shortening the window might appear to be a solution to this problem, but this leads to a serious degradation of the coherency value for weak reflections where the signal-to-noise ratio is low.
If high resolution velocity estimates are desired, they may be obtained by computing the velocity spectra on a fine grid of velocities and times. This increases the computational costs greatly. It is therefore desirable to have a method for automatic velocity picking that has a higher resolution without the computational costs that would be incurred in computing the velocity spectra on a fine grid of velocities and times. One method that has been used in the prior art starts out with a computation of the velocity spectra on a coarse grid. This is followed by a search technique, using steps determined from a Fibonacci series, to compute the coherency at additional trial velocities. The Fibonacci search technique would be familiar to those skilled in the art. By successively reducing the step size, the velocity that gives a maximum of the coherency for a time window is determined. The coherency is calculated only at the velocities determined by the Fibonacci series rather than on a fine grid of velocities between the bounds specified by the analyst. This can lead to a reduction in computation time.
The Fibonacci search method still has the smearing caused by the finite time window, and the effects of stretching are present. In addition, this search technique is strictly one-dimensional, i.e., the zero-offset times are treated independently of each other and without regard to whether in fact there is a reflection within the window.
From the foregoing, it can be seen that an improved method for automatic picking of NMO velocity functions is needed. Such a method should provide for automatic picking of the NMO velocity function for traces within a CMP gather within bounds specified by the analyst. It should also provide a high resolution estimate of the NMO velocity function without an excessive computational burden, and should provide this high resolution estimate without the stretching and smearing introduced by computing the coherency function over a finite time window. The present invention satisfies these needs.